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Electron energy level statistics in graphene quantum dots 
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Motivated by recent experimental observations of size quantization of electron energy levels in graphene 
quantum dots [7] we investigate the level statistics in the simplest tight-binding model for different dot shapes 
by computer simulation. The results are in a reasonable agreement with the experiment which confirms 
qualitatively interpretation of observed level statistics in terms of "Dirac billiards" without taking into ac- 
count many-body effects. It is shown that edge effects are in general sufficient to produce the observed level 
distribution and that even strong bulk disorder does not change the results drastically. 
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Graphene attracts enormous attention now, due to 
its interest both for fundamental physics, such as oppor- 
tunities to simulate in condensed matter experiments 
subtle quantum relativistic effects and for potential ap- 
plications, as a planar, high-mobility material for "post- 
silicon" electronics [TJ EJ 02 Efl [5] . Graphene-based nan- 
odevices are subjects of especial interest. Recently, size 
quantization effects have been observed in graphene 
nanoribbons [6] and quantum dots (QD) [7]. It was 
demonstrated that for QD smaller than 100 nm the 
electron energy spectrum is essentially irregular demon- 
strating a "chaotic" behavior. The latter can be dis- 
cussed in terms of random matrix theory for a single- 
particle problem [8l[9l[T0]. However, in general, because 
of interplay of size quantization and Coulomb blockade 
correlation effects can be important for QD pj]; more- 
over, in some limiting cases the "chaotic" energy spec- 
trum can be described even purely classically, in terms 
of only Coulomb energies [12]. Therefore, to understand 
properly the experimental data [7j some theoretical ef- 
forts are necessary. 

In this Letter we present the results of straight- 
forward computer simulations of level statistics in 
graphene QD, adopting the simplest onc-clcctron pic- 
ture, up to 16000 sites (which is comparable to the sizes 
of smallest QD investigated experimentally). It turns 
out that already this approach does allow to reproduce, 
in a semi-quantitative way, the observed energy distri- 
bution. 

As the model Hamiltonian, we take the simplest 
nearest-neighbor tight-binding model on a hexagonal 



120 
100 
80 



& 60 



40 
20 




x ) e-mail: M.Katsnclson@science.ru.nl 



3 
S 



Fig. 1 (Color online) Level-spacing distribution P(S) 
as obtained from experiment [7] for the case of 40 nm 
QD (The raw experimental data are by courtesy of K. 
Novoselov and A. Geim). 



lattice as introduced in Ref. [13], enclosed in some ge- 
ometrical shape, such as a circle, triangle, etc. The 
Hamiltonian reads 

F = -*E E (4 A* + b J,* a i,*) 
(i,j> °-=T4 

+E E l • E E l AA^ 
+E E v i a L%a+J2 E v i b hh*> w 

i cr=t,J. j <T=t,l 

where t is the nearest-neighbor hopping energy, aj 
(a t a ) creates (annihilates) an electron with spin a on 
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one of the sub-lattices of the hexagonal lattice and bj a 
{bj a ) creates (annihilates) an electron with spin a on 
the other sub-lattice of the hexagonal lattice. In the 
expression for the hopping term, the indices i and j run 
over all nearest neighbors only. 

Following Ref. [T3], we also consider the effects of a 
staggered on-site potentials Vi and Vj that alternate as 
we move along the boundary B of the lattice; as shown 
in Ref. [14] , the continuum limit with a very large Vi and 
Vj corresponds to the "infinite mass" boundary condi- 
tion [8]. The physical origin of this potential can be 
related to, for instance, the magnetic moments at the 
zigzag edges [TSJ [16] . Of course, in a simulation we can 
easily study the effect of Vi by considering the cases 
V = Vj = and V t ^ and Vj ^ 0. Optionally, 
to study the effect of disorder, we use uniform pseudo- 
random numbers to choose the hopping integrals from 
the interval [t — 5,t + 6] and/or we add on-site potentials 



and Vj in the range [- 
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The geometrical shape "cuts out" a piece of the 
hexagonal lattice such that armchair and zigzag parts of 
the boundary appear (for a general discussion of bound- 
ary conditions for the tight-binding model of graphene, 
see Ref. [H]). Furthermore, we study the effect of disor- 
der that results from removing hexagons from the reg- 
ular hexagonal layer. 

The eigenvalues of Hamiltonian Eq. (TT]) are obtained 
by exact numerical diagonalization. As we are inter- 
ested in the part of the spectrum that, in the continuum 
limit corresponds to the spectrum of the Dirac equa- 
tion, we limit the search for eigenvalues to the interval 
[— 0.4i, OAt] . In our numerical work, we adopt units such 
that H = 1 and we set the hopping integral t = 1. 

For reference, in Fig.[TJ we show the experimental re- 
sults of the level spacing distribution for a graphene dot 
of 40 nm diameter [7] . The dimensionless level spacing 
S is defined as the energy difference AEi = Ei — Ei—i 
between successive levels, divided by the average (AEi) 
of the energy differences between successive levels. The 
number P(S) gives the number of energy difference for 
which S - A/2 < AEi/{AEi) < S + A/2, where A is 
the bin size of the histogram. The experimental data for 
the distribution have been shifted by 0.02 V to remove 
an ambiguity in the definition of S = 0. 

For comparison, Fig. [2] shows the level distribution 
for a hexagonal lattice with periodic boundary condi- 
tions, a lattice that has not been cut-out using some 
geometrical shape. It is clear that its spectrum does 
not resemble the one observed experimentally. 

Next, we consider the case of a circular dot. In the 
continuum approximation with the infinite mass bound- 
ary conditions [8] this case is special, with separable 




Fig. 2 (Color online) Level-spacing distribution P(S) for 
a hexagonal lattice with periodic boundary conditions, 
using 170 different values for both k x and k y . 



variables and a regular energy spectrum. In reality, al- 
though we use a circle to cut out from the infinite lattice, 
the boundary of this lattice is irregular: There are short 
and long pieces of the armchair and zigzag boundaries. 
This kind of irregularities destroy the circular symmetry 
of the the geometrical shape completely. 
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Fig. 3 (Color online) Level-spacing distribution P(S) for 
a hexagonal lattice bounded by a circle. The number 
of hexagons inside the circle is 8202. 



In Figs. [3] and HI we depict the results for P(S) and 
the density of states DOS(£') for the case of the perfect 
hexagonal lattice, bounded by a circle. Although there 
is a clear background linear dependence of the DOS(E') 
on E, it is also clear that there are fluctuations due to 
size quantization. In the DOS, the peak around zero en- 
ergy is due to the existence of the zigzag edge states, as 
expected for a generic boundary |14j . Comparing Fig. [3] 
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Fig. 4 (Color online) Density of states DOS(i5) as a 
function of the energy E for a hexagonal lattice con- 
fined to a circle. The number of hexagons inside the 
circle is 8202. 
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Fig. 5 (Color online) Same as Fig. [3] except that the 
staggered potential at the boundary edges V = 100, 
the hopping integrals fluctuate by maximum 20% (S — 
0.2) and the on-site potentials fluctuate in the range 
[—0.2, 0.2] The number of hexagons inside the circle is 
8484. 

with the experimental result Fig. [TJ we conclude that 
there is little resemblance. 

In contrast, by introducing various forms of disorder 
to the same system, we find semi-quantitative agree- 
ment, as shown in Figs. [5] and [H Our numerical ex- 
periments (not all results shown) suggest that the pres- 
ence of an alternating boundary potential can change 
the qualitative features of P(S) significantly. Including 
various forms of disorder (<5 ^ 0, v ^ and remove 
some hexagons), the level distribution P(S) (see Fig. [6]) 
looks similar to the experimental result (see Fig. [T]). Of 
course, we cannot expect quantitative agreement: The 
number of items in the experimental data is about 60 



Fig. 6 (Color online) Same as Fig. [3] except that the 
potential at the boundary edges V = 10, about 20% of 
the hexagons have been removed, the hopping integrals 
fluctuate by maximum 20% (S — 0.2) and the on-site 
potentials fluctuate in the range [—0.2, 0.2] The number 
of hexagons inside the circle is 6540. 



while in the numerical simulations there are about 600 
eigenvalues that contribute to P{S). This may explain 
why the simulation results for P(S) show a more ex- 
tended tail than the experimental result for P(S). 

Finally, to study the effect of the shape of the bound- 
ary on the level spacing distribution, we have calculates 
P(S) for various lemon-shaped billiards (results not 
shown) defined by y — ±(1 — \x\ d ) for x G [—1, 1] [T7]. As 
a function of the shape parameters d > 0, these billiards 
look like a square (d = 1, oo) or two parabola {d — 2) or 
have some intermediate lemon-like shape. It has been 
shown that, as a function of d, the these billiard, clas- 
sical as well as quantum mechanically, exhibit regular 
and chaotic behavior [17] • If Figs. [7] and we present 
some typical results for a quarter lemon with d = 3.1. 
Our motivation for presenting the results of a quarter 
lemon instead of the complete lemon is to show a case 
in which there are long stretches of armchair and zigzag 
boundaries (the edges along the x and y axis) and an 
irregular boundary (the curved edge). Although on pur- 
pose, we did not include the alternating potential at the 
boundary sites, the disorder resulting from the irregular 
shape together with the fluctuating on-site potential v 
seem to be sufficient to observe a P(S) that is similar 
to the distribution observed experimentally. Compar- 
ing Figs. 7 and 8 one can see that the level statistics 
is not too sensitive to the bulk disorder. There is only 
randomness due to boundaries themselves in one case 
(Fig. 7) and a random potential v is introduced, addi- 
tionally, in the other one (Fig. 8), but the results look 
very similar. 
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Fig. 7 (Color online) Level-spacing distribution P(S) 
as a function of the level spacing S for a hexagonal lat- 
tice bounded by a quarter of a lemon-shaped billiard 
for V = S = v = 0. Parameter of the lemon: d = 3.1 
(d — 1 is half of a triangle, d = 2 is a half of parabola). 
The number of hexagons inside the billiard is 7861. 
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Fig. 8 (Color online) Same as Fig. [7] except that v is 
chosen at random from [—0.2,0.21. 



To conclude, it seems that disorder due to random- 
ness of the edges is, in principle, enough to reproduce 
the experimentally observed level distribution which 
makes the term "chaotic Dirac billiard" quite reason- 
able. One may need the local on-site disorder and some 
disorder in the hopping integrals to get semi-qualitative 
agreement with experiment but, on the other hand, the 
simulated systems do not have the same shape and are 
not as large as the experimental ones so a complete 
quantitative agreement is, anyway, hard to expect. 

It is worth mentioning that the continuum approxi- 
mation may be a bit dangerous when discussing the level 
statistics in graphene QDs. In our simulations, we do 
not see any essential differences between "regular" (cir- 



cular) and "irregular" (lemon-shaped) billiards. Even 
for tens of thousand sites the edge is essentially irregular 
with staggered armchair and zigzag pieces. Therefore, 
even for circular quantum dot we may have a "chaotic" 
energy level distribution. 

We are thankful to Andre Geim and Kostya 
Novoselov for helpful discussions and for providing us 
with the original experimental data used in Fig. [1] and 
to Boris Shklovskii for illuminating discussions. The 
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Onderzoek der Materie (FOM) (the Netherlands) 
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